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Abstract 

We consider a threshold-crossing spiking process as a simple model 
for the activity within a population of neurons. Assuming that these 
neurons are driven by a common fluctuating input with Gaussian 
statistics, we evaluate the cross-correlation of spike trains in pairs 
of model neurons with different thresholds. This correlation function 
tends to be asymmetric in time, indicating a preference for the neu- 
ron with the lower threshold to fire before the one with the higher 
threshold, even if their inputs are identical. The relationship between 
these results and spike statistics in other models of neural activity are 
explored. In particular, we compare our model with an integrate-and- 
fire model in which the membrane voltage resets following each spike. 
The qualitative properties of spike cross-correlations, emerging from 
the threshold-crossing model, are similar to those of bursting events 
in the integrate-and-fire model. This is particularly true for general- 
ized integrate-and-fire models in which spikes tend to occur in bursts 
as observed, for example, in retinal ganglion cells driven by a rapidly 
fluctuating visual stimulus. The threshold crossing model thus pro- 
vides a simple, analytically tractable description of event onsets in 
these neurons. 



1 Introduction 

Probing the relationship between a stimulus and the spike train, generated 
by a population of neurons, is a central theme in the study of neural activity 
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(Rieke, Warland, Stevennick, & Bialek, 1996). Correlations in the spike 
timing of different neurons are of interest, in this context, for two main 
reasons: First, spike correlations are informative about the structure of the 
neural code, beyond what could be inferred from the firing properties of 
single neurons alone. Second, correlations are often thought to reflect the 
structural properties of the neural network (Ginzburg & Sompolinsky, 1994). 

Correlated firing in a pair of neurons can arise because of two distinct 
reasons: the existence of a connection between the neurons (direct or indi- 
rect), and the existence of a common input to the two neurons that is varying 
in time. Here we focus on the second possible source for correlated activity. 
We consider model neurons that receive an analog, continuous stimulus, and 
respond to it by generating a discrete sequence of spiking events. The main 
question that we address is how the statistical properties of the fluctuating 
stimulus affect the structure of spike correlation functions. 

Correlated inputs to neurons were considered theoretically mainly in 
context of their effect on single-neuron properties, such as the firing rate 
(Tuckwefl, 1988b; Salinas & Sejnowski, 2000; Moreno, Rocha, Renart, & 
Parga, 2002; Kuhn, Aertsen, & Rotter, 2003), the coefficient of variation 
(CV) (Tuckwell, 1988b; Brunei & Sergi, 1998; Salinas & Sejnowski, 2000; 
Stroeve & Gielen, 2001; Salinas & Sejnowski, 2002; Schwalger & Schimansky- 
Geier, 2008), and the spike-triggered average stimulus (Kanev, Wenning, & 
Obermayer, 2003; L., Gerstner, Sz Richardson, 2006; Paninski, 2006). These 
works considered integrate-and-fire model neurons (Tuckwell, 1988a; Gerst- 
ner & Kistler, 2002). 

The non-leaky integrate-and-fire model is relatively tractable analytically 
(Tuckwell, 1988b; Paninski, 2006). In comparison, analytical treatment 
of the leaky integrate-and-fire (LIF) neuron is considerably more difficult 
(Burkitt, 2006). Some analytical results are available for the firing rate 
(Salinas & Sejnowski, 2000; Moreno et al., 2002; Kuhn et al., 2003), whereas 
analytical results for the ISI and the CV are available only in particular 
limits such as slowly varying (Moreno-Bote & Parga, 2006; Schwalger & 
Schimansky-Geier, 2008), or binary (Salinas & Sejnowski, 2002), inputs. 
Other results were obtained for the LIF model from computer simulations 
(Salinas & Sejnowski, 2000; Stroeve & Gielen, 2001). In particular, Stroeve 
and Gielen (2001) included a simulation study of correlations in spiking 
of LIF neurons that receive partly overlapping input. Correlations due to 
partly overlapping input were also studied analytically, but only in the limit 
where the input is slowly fluctuating in time (Moreno-Bote & Parga, 2006). 

Here we consider a simpler model of neural response, where neurons 
spike whenever a generating potential g(t), linearly related to the neuron's 
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stimulus, crosses a threshold in its rising phase. For sufficiently simple 
stimuli this model is analytically tractable, which allows for spike correlation 
functions to be evaluated in closed form. 

Threshold crossing processes without a reset were previously analyzed 
as models for neural firing. The spike auto-correlation function of a neuron 
was considered by Jung (1994), while making specific assumptions on the 
nature of the fluctuating potential. Here we evaluate the auto- and cross- 
correlation in spike timing of neurons with different thresholds, while making 
fewer assumptions on the generating potentials eliciting the spikes. These 
are assumed to be Gaussian, and may be identical or partially overlapping 
in different neurons. 

As in the case of Linear-Nonlinear (LN) models (Korenberg & Hunter, 
1986; E. Chichilnisky, 2001), we assume that a neuron responds to a tem- 
poral convolution of its stimulus with a linear kernel. However, in the LN 
model spiking is stochastic, whereas in the threshold-crossing model the 
spike timing is precisely determined by the stimulus. This aspect of the 
model is motivated by the observation, in various neural assays, of responses 
that are precisely repeatable across multiple trials (Mainen & Sejnowski, 
1995; Berry, Warland, & M., 1997; R. R. de Ruyter van Steveninck, Lewen, 
Strong, Koberle, & Bialek, 1997; Meister & Berry, 1999; Keat, Reinagel, 
Reid, & Meister, 2001; Uzzell & Chichilnisky, 2004), much more than can 
be described by Poisson statistics, particularly when the stimulus is strongly 
fluctuating in time. 

The model and our main assumptions are presented in Sec. II. Before 
discussing spike correlation functions, we first evaluate the firing rate (Sec. 
Ill), the spike-triggered average stimulus, and the spike-triggered covariance 
(Sec. IV). We then evaluate spike auto- and cross-correlation functions 
(Sec. V). These results are compared with computer simulations of model 
neurons that spike according to several variations of the leaky integrate-and- 
fire model (Sec. VI). 

2 Formalism and assumptions 

We consider one or more neurons that respond to the same generating po- 
tential g and elicit a spike whenever g crosses a threshold in its rising phase. 
The spike train generated by the neuron i can be written as 

Xi(t) = 5[g(t)-e i }g(t)G[g(t)} (1) 
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where g is the generating potential. The Heaviside step function 6 [g(t)] 
restricts the firing to the upward crossing events. 

We assume that g is stationary, Gaussian, and has zero mean. The 
generating potential is thus fully characterized by its correlation function, 

(g(t)g(t')) = w(\t'-t\) (2) 

where the brackets () stand for an ensemble average over all possible real- 
izations of the fluctuating generating potential g. Later on, in Sec. V, we 
consider a more general situation of a population of neurons, whose gener- 
ating potentials gi are jointly Gaussian, characterized by their covariance 
functions 

(g i (t)g j (t')) = w ij (t'-t). (3) 

The properties of a model neuron, responding to a stimulus whose mean dif- 
fers from zero, can be obtained from the zero-mean case simply by adjusting 
the threshold. 

For a single neuron, the behavior of w(At) at small At determines the 
firing rate and the short-time behavior of the spike auto-correlation (Sec. V). 
We assume that an expansion of w(At) exists around At = 0, 

w(At) = W + W!\At\ + ^W 2 At 2 + ^VF 3 |At| 3 + . . . (4) 

We further assume that W± vanishes and that Wi is negative, for reasons that 
will become clear later on. As seen below, it is necessary to treat separately 
different classes of processes, depending on which other coefficients in the 
expansion are non-zero. One important class is the case where W3 7^ 0. This 
is the typical situation if a causal filter is involved in the generation of g(t). 
For comparison we also briefly consider, in Sec. V, the case where all the 
coefficients with odd indices vanish; for example, a process with a Gaussian 
correlation function, w(t) = <r 2 exp(— t 2 /2r 2 ). 

Relationship between g and the stimulus 

The generating potential in our model is linearly related to the stimulus s, 

g = fo S (5) 

where the symbol 'o' stands for convolution. We assume that s is a zero- 
mean, uncorrelated Gaussian process: 

(s(t)s('t)) = a 2 S(t'-t) (6) 
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The correlation function of g is then 

w(At) = a 2 J f(t)f(t + At)dt (7) 

Note that Wo = a 2 , J™ [f(t)] 2 dt > and that if / is sufficiently regular, 

W 2 = -aifZ o lf(t)} 2 dt<0. 

The assumption that s is uncorrelated is made for simplicity of the pre- 
sentation, and because uncorrelated flickering stimuli are often used exper- 
imentally The results below generalize in a straightforward manner also to 
the case where s is a correlated Gaussian process. 



Specific forms for / and w 

In the following sections we first derive results that are valid generally, for 
any w(At). We then illustrate these results with a specific example, where 
we assume a particular form of /. In these examples the filter / is a causal 
filter of the form 

'0 t < 

f(t)=) e -*M_ e -tAi ^ Q ( 8 ) 

T2 - n 

which is a combination of two single-exponential filters with time constants 
Ti and T2- Using Eq. (7) 

w{At) = fJ 2 r 2exp(-|At|/r 2 ) -nexp(-|At|/ri) ^ 

^2 - n 



where 



a 2 = - a ° (10) 
2(ri+r 2 ) V ' 

is the variance of g. 

In the particular case where n = ti = r, / is an 'alpha' filter, 

r o t < o 

f{t) = { ^ t > o < u » 

and 

2 

w(At) = — (\At\ +r)e- |A * l/r (12) 
r 

The expansion of w(At) for small At, Eq. (4), yields W\ = and: 



W ^ W ^ W Tl + T2 - 2 W T l + TlT2 + T 2 ^ 
Wo = a , W 2 = , W 3 = — n-^-cr , W4 = ^ cr 

T1T2 rfr| rfr| 



(13) 
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3 Firing rate 



We begin with the relatively simple problem of evaluating the firing rate of 
a single neuron (Rice, 1954), 

r = ( X (t)) (14) 

where \ is given by Eq. (1). Although this quantity has been calculated 
before, we derive it here in some detail, because the derivation general- 
izes to the higher-order moments that are calculated later (whose detailed 
derivation is presented in the appendices.) The firing rate r can be written 
as 

/■oo 

dqq-p(9,q) (15) 



where p(g, q) is the joint probability distribution for the generating potential 
to be equal to g and for its derivative, at the same time, to be equal to q. 
To evaluate this and similar quantities, we use an identity that holds for a 
general Gaussian signal g{t) with correlation function w(At): If £i . . . £ n are 
all scalar random variables that depend linearly on the process g, 



= J dtOi(t)g(t), 



then the joint probability distribution of £i, • ■ ■ , C, n is equal to 



p(Ci) • • • >Cn) = ^exp 



2 



(16) 



(17) 



where Z = [(27r)™detA] 1/2 and 

^ = (18, 

The generating potential and its derivative at time t = correspond to 
Ci = 5(0) and C2 = 5(0), therefore a\ = 5{t) and a 2 = —S(t), so that 



A 



W 






-w 2 



(19) 



The off-diagonal terms are proportional to u/(0), which vanishes if W\ = 0, 
as we assumed 1 . We thus have 

„2 „2 



1 



; exp 
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2TT -2W2) 



(20) 



1 li Wi / 0, w(At) does not have a derivative at At = 0, and the derivation performed 
here fails. A typical example is given by uncorrelated noise passed through a single- 
exponential filter. For this process the firing rate diverges (see below) . 
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and the integral (15) yields 



2tt [ W \ 6XP V 2W ) 



1 \-w 2 ] 1/2 f 9 2 \ 



(21) 



Hence W2 must be negative for the firing rate to be finite. 
For the generating potential of Eqs. (8) and (9), 



r "2vr(r 1 r 2 ) 1/2eXP ( ' 



(22) 



If either t\ or T2 vanish, the firing rate diverges. Hence a generating potential 
obtained by convolving uncorrelated noise with a single exponential filter 
has an infinite rate of threshold crossings. In such a process there are finite 
intervals in which the number of crossings is infinite, and others in which no 
firing occurs 2 . 

4 Spike triggered average stimulus 

The spike triggered averaged stimulus (STA) is the mean stimulus given that 
a spike was generated at a particular time. Because the input is assumed to 
be stationary, the STA is only a function of the time difference relative to 
the spike time, 



where r is the firing rate, x 1S determined by g, as in Eq. (1), and the 
spike time was arbitrarily chosen to be zero. If one assumes that spiking is 
described by a LN model with filter /, and if the stimulus is uncorrelated 
Gaussian noise, /sta(A£) is proportional to /(—At) (E. Chichilnisky, 2001; 
Paninski, 2003). Hence it is common to probe the spatio-temporal filter 
applied by a neuron on its stimulus by measuring /sta (Schwartz, 2006), in 
particular in the visual sensory system (Rieke, 2001; E. J. Chichilnisky & 
Kalmar, 2002; Baccus & Meister, 2002; Rust, Scwartz, Movshon, Sz Simon- 
celli, 2005; Hosoya, Baccus, & Meister, 2005). 

To compute /sta for a threshold- crossing spiking neuron we need the 
joint probability distribution function p(gi, q±; S2) for g(0) = <?i, g(0) = qi, 
and s(At) = S2, evaluated in Appendix A. The STA is then given by 



This is related to the short-time properties of the spike auto-correlation function 
(Sec. V) in the limit where n or T2 — » 0. 



WAt) = - (s(At) ■ X (0)> 



(23) 




(24) 
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where A = W /a 2 = f 2 {t)dt and B = -W 2 /a 2 = f 2 (t)dt. 



The first term in Eq. (24) is proportional to the expectation for the STA 
in a rate-based linear-nonlinear model (E. Chichilnisky, 2001). However, for 
a threshold-crossing spiking neuron it correctly describes the STA only in 
the limit of large 9: in contrast, when = 0, the STA is proportional to 



The result that /sta is a linear combination of the filter and its derivative 
has a simple geometric interpretation, which is most easily seen by thinking 
of the stimulus history as an n-dimensional discrete vector. The condition 
for spiking at time t = involves g and its derivative, which are the inner 
products of the stimulus history with two vectors, / and /'. Because s is 
uncorrelated, its history vector is distributed in a radially symmetric manner 
in n-dimensional space. The STA must lie in the two-dimensional sub-space 
spanned by the vectors / and /'. 

This result is illustrated in Fig. la, for the generating potential of Eqs. (11) 
and (12). With the threshold varying from to 2.5a, the shape of /sta varies 
significantly. The dotted lines in the figure (overlapping with the solid line) 
were obtained from a discrete simulation, using a time step dt = r/100 
and averaging over a simulation run T/r = 10 6 for 9 = and 6 = 1, and 
T/t = 10 8 for 6 = 2.5. 

It is interesting to compare these results with what is expected from a 
leaky integrate-and-fire model neuron (Tuckwell, 1988a; Gerstner & Kistler, 
2002). For this case there is no known analytic form for the STA (see, 
however, Kanev et al., 2003; L. et al., 2006; Paninski, 2006), but we can 
evaluate the STA numerically. We assume that the membrane potential of 
the model neuron is related to the input current via 



with a reset of the membrane potential to u r whenever u reaches the thresh- 
old 9, whereas the input current is related to the stimulus s by 



The sub-threshold dynamics of u depends on the stimulus in the same way 
as the generating potential g depends on stimulus in the threshold-crossing 
model. The difference between the two models lies in the existence of a 
reset. For simplicity, we set u r = and t\ = t 2 for the rest of this section 
(see also Sec. 6 and Appendix D). 



f(At). 




du 



(25) 




(26) 



8 



Figure lb shows a comparison between /sta in the integrate-and-fire 
model (simulation, dotted line: dt/r = 0.01 and T/r = 10 6 ) and /sta in the 
threshold- crossing spiking model [Eq. (24), solid line.] The parameters are 
as in panel a, with a threshold 9 = 2.5a. The two models agree very well 
in their prediction and differ, significantly, from that of a LN model, where 
/sta (At) oc f(-At) (dashed line). 

For a lower threshold, = 0.25, the threshold-crossing spiking model 
and the integrate-and-fire model no longer yield similar predictions for /sta 
(solid and dashed lines, Fig. lc). The difference between the cases of high 
and low threshold is possibly related to the typical inter-spike interval which, 
in the case of a high threshold, is longer, allowing the input to decorrelate 
between subsequent spikes. The statistics of the input before a spike are 
thus less influenced, in the case of high threshold, by whether a reset has 
occurred following the previous spike. 

Finally, we note that the threshold-crossing model is similar to a two- 
dimensional LN model involving two filters, equal to /(— t) and /'(— t). In 
the analogous LN model the nonlinear transfer function must be chosen to 
be narrowly peaked at values of the generating potential close to the thresh- 
old, and an appropriate dependence on the derivative must be included as 
well 3 . The result that the STA is a linear combination of the filter / and 
its derivative thus agrees with the general property of multi-dimensional LN 
models, that the STA lies within the subspace spanned by the linear kernels 
determining the spiking rate (Paninski, 2003; Schwartz, 2006). Similarly, we 
may expect the spike-triggered covariance (R. de Ruyter van Steveninck & 
Bialek, 1988; Schwartz, 2006) to reveal the linear filter / and its derivative 
/' as spanning this sub-space. This is indeed the case (Appendix A.) 

5 Spike correlations 

We next consider correlations in the spike timing of two neurons, assuming 
that they receive identical inputs, but possibly differ in their thresholds 
0± t 2- Because the generating potential is assumed to be stationary, the spike 
correlation function depends only on the time difference between spikes: 

c(At) = (xi(0)-X2(At)> (27) 

From the definition of xi,2, Eq. (1), the quantity inside the brackets depends 
on the generating potential and its derivative at t = and at t = At. 

There is a difference between the two models, however, in that the two-dimensional 
LN model produces a variable number of spikes within the narrow spiking event. 
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The detailed form of c(At) is derived in appendix B, and here we present 
the main results. The spike correlation function can be written as 



/»oo /»oo 

c(At) = p(9 1 ;9 2 ) ■ Ml dq 2 qi-q 2 
Jo Jo 



(28) 



(detm- 1 ) 172 



2tt 



exp 



i (q - vof m 1 (q-qo) 



where p(9\;9 2 ) is the joint probability distribution for #(0) = 9\ and g(At) = 
#2 [see Eq. (66)], and q T = (qi,q 2 )- The quantity in the second line of the 
equation is the conditional probability distribution function for </(0) = q\ 
and g'(At) = q 2 , given that g(0) = 9\ and g(At) = 9 2 . This conditional dis- 
tribution is Gaussian and is characterized by its mean qo and the covariance 
matrix m, whose values are derived in the appendix [Eqs. (64) and (65)]. 



The spike auto-correlation function 

We briefly discuss the spike auto-correlation function, corresponding to the 
case 9\ = 6 2 = 9, and focusing on the behavior at small At. Because the 
spike train x(t) is a point process, its auto-correlation function necessarily 
includes a contribution 



c(At) = r5{At) + ... (29) 

The discussion here concerns nonzero values of At, i.e., the occurrence of 
spikes in addition to the one at t = 0. The behavior for small At is remark- 
ably different for the two classes of generating-potential statistics discussed 
in Sec. II: If W3 7^ 0, as in the example of Eq. (13), the spike auto-correlation 
function tends to a finite value when At — > 0: 

W 3 ( 9 2 \ 

G{At) A^-SWoWtf/** 3 * I." 2W~ ) (30) 

It is interesting to look at this quantity divided by r, which represents the 
firing rate at time t = At after the occurrence of a spike at t = 0: 

cm _ w. (31) 

r -2wV3W 2 

This ratio does not depend on 9. By comparing with r itself we see that, 
if the threshold is high, the spike-conditioned firing rate is much higher 
than the average firing rate. In other words, once the neuron has fired, it 
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is likely to soon fire again, compared to its baseline firing rate. When the 
threshold is zero, the spike-conditioned firing rate can be larger or smaller 
than r, depending on the expansion coefficients Wi and W3. This result is 
illustrated in Fig. 2a, where w is given by Eq. (12) 4 . 

We note that when W3 / as above, the second derivative of the gener- 
ating potential is an unbounded random process, and the first derivative of 
g is not smooth. This facilitates the generation of successive spikes at short 
intervals, since the derivative of g is required to be positive at the interval's 
edges, and negative somewhere in between. 

In contrast, when the first derivative of g is a smooth process, we may 
expect the occurrence of spikes at vanishingly small intervals to be consid- 
erably less likely. Indeed, if w(At) has no irregularities at At = 0, so that 
all Wi with odd indices vanish, the occurrence of spikes separated by short 
intervals is strongly suppressed, c(At) ~ (At) 4 [Appendix C, Eqs. (84) and 
(87)] . This is illustrated in Fig. 2b for a generating potential with a Gaus- 
sian correlation function, w(At) = cr 2 exp(— At 2 /2r 2 ). Note that Wo and 
W2 are the same in the two parts of Fig. 2; consequently, the mean firing 
rates are identical in these two cases. 

Spike cross-correlation 

We next consider the spike cross-correlation function of two neurons with 
different thresholds, firing in response to the same stimulus. We note, first, 
that unless 9\ = 62 the correlation function is not symmetric with respect 
to replacement of At by —At or, equivalently, with respect to exchange of 
the two neurons {9\ <-> 62)- The generating potential is required to have a 
positive derivative at t = and at t = At, and therefore we may expect a 
higher probability for joint spiking if the higher of the two thresholds is set 
at the later time (see Fig. 3) . This is indeed the case - as demonstrated in 
Fig. 4a for the generating potential described by Eq. (12), and with thresh- 
olds set as 9\ = O.80" and 62 = o. The preference for one neuron to fire after 
the other neuron is particularly prominent at short time scales of order r, 
but it also has a signature at large time scales (Appendix B). 

It is instructive to compare the threshold-crossing model with the pre- 
diction of a one-dimensional linear-nonlinear (LN) model (E. Chichilnisky, 
2001). For this comparison, we assume that the linear filter of the LN model 
is the same as the filter / in the threshold-crossing model, and that in both 

4 Note that the conditional firing rate diverges if n — > or if T2 — » 0, which can be 
seen from Eq. (13). This is consistent with the divergence of the average firing rate in this 
limit, and the existence of finite time intervals in which the number of spikes is infinite. 
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models the stimulus is uncorrelated Gaussian noise. In the LN model the 
spike correlation function depends also on the choice of the non-linear func- 
tion applied to the outcome of the linear filter. Two particular choices are 
made in the examples shown in Fig. 4b: The first is linear rectification, 
4>{x) oc Q(x — 9) where 6(x) = x for x > and Q(x) = for x < 0, with 
thresholds 0\,0 2 chosen as in Fig. la (solid line). 

The second example (dashed line) is the case where <f)(x) is non-zero only 
in a narrow range around 6, 4>{x) oc S(x — 0). Mathematically, this example 
compares more directly with the threshold-crossing model, because the LN 
model produces spikes only at the threshold crossings. The correspondence 
can be seen in the spike cross-correlation function c(At) obtained from the 
LN model in this case, which is calculated in Appendix E [Eq. (100).] In 
the limit of large At, c(At)/rir 2 — 1 in the LN model [Eq. (104)] is equal 
to the first term in the large At expansion for the threshold-crossing model, 
Eq. (74). 

We note two important qualitative differences between the LN model 
and the threshold-crossing spiking model: First, in the LN model, the spike 
correlation function is symmetric under time reversal: c(— At) = c(At) 5 . A 
second qualitative difference between the two models, seen in Fig. 4, is that 
the firing of the two neurons is less correlated in the LN model, compared 
to the threshold crossing spiking model, at short time scales. 



The limit At -> 

In the limit of small At the behavior of c(At) is substantially different if the 
neuron with the higher threshold fires before or after the neuron with the 
lower threshold 6 . To simplify notation we assume here that 2 > 9\. 

At > 0: For small At the spike correlation function scales with At as 



c(At) 



(9 2 - 0i 
At 3 



-exp 



1 



2W 2 



At 



+ 



7i_ 
At 



(32) 



where 



7i 



W 3 
"6W 2 2 



(02 ~ 0lf 



(33) 



5 This property arises for the following reason. The spike cross-correlation at any two 
times t\ , t2 depends on the joint probability distribution of the generating potential values 
at these times. Assuming that the generating potential is Gaussian, this joint distribution 
function is necessarily symmetric. 

6 A weak asymmetry between positive and negative At exists also in the limit At — » co. 
This limit is considered in Appendix B. 
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This approximation, with the appropriate prefactors taken from Eqs. (60), 
(76), and (78) is plotted in Fig. 4a (dashed line) for the particular input 
considered in this figure. 

As At — > 0, the exponential decay in Eq. (32) wins over the ~ (At)- 3 di- 
vergence for sufficiently small At. At larger values of At the spike correlation 
function peaks, and then decays algebraically. Hence the neuron with the 
larger threshold tends to fire after the neuron with the lower threshold, with 
a typical latency given by the position of the peak. If we ignore the term 
71 /At, which is legitimate when the difference in thresholds is sufficiently 
small, the maximum is at At = At* where 

Af = TW (34) 

For the generating potential of Eqs. (8) and (9), 

A t - . (35) 

Equation (34) indicates that the most likely latency increases linearly 
with 62 — 9\. An illustration of this result is shown in Fig. 5, where / and 
w are as defined by Eqs. (11) and (12), 62 = 0.5a, and c(At) is plotted for 
three values of 0\: 0.2a, —0.2a, and — 0.5a (solid lines). With increase in 
82 — 9\ the peak of the spike correlation function becomes wider, and occurs 
at larger latencies. The arrows indicate the prediction for the position of 
the peak, Eq. (35), which matches the actual position very well even when 
62 — 61 is relatively large. 

At < 0: In this case the process g is required to have a positive derivative 
at t = At < 0, a negative derivative within the interval (At, 0) and, again, a 
positive derivative at t = 0. Such a trajectory becomes increasingly unlikely 
when At — > 0, and the spike correlation function decays to zero. The leading 
contribution is of the form 

<«> 

if W3 / 0. If W3 = W5 = the decay is even stronger, logc(At) oc (At) -6 
(Appendix C.) Equation (36) introduces a characteristic time scale for the 
inhibition, scaling as 



hf /3 ( rlrl N 1/3 -A~\ 2 / 3 



(37) 



where the expression on the right hand side holds if w(At) is given by Eq. (9). 
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Partly overlapping inputs 

We next generalize the analysis to the more realistic case where the inputs 
to different neurons are only partially overlapping: 

ft = / o (s + &) (38) 

The noise inputs £j are assumed to be zero-mean, jointly Gaussian processes 
that are uncorrelated with s and with each other: 

(UtMt')) = aa 5 ij S(t' -t). (39) 

Here a is the ratio between the standard deviation of the noise and that of 
the common stimulus. The covariance of the generating potentials is then 

Wij(At) = w(At)(l + a 2 <%) (40) 

The spike auto- and cross-correlation functions can be evaluated in a similar 
manner as for the noise-free case (Appendix B). 

As an example, we consider the 'alpha' filter, Eq. (11). In Fig. 6 the 
standard deviation a of the common stimulus is kept fixed, while a is var- 
ied. Figure 6a shows the spike correlation function of two model neurons 
with the same threshold 9 = a and with weak independent noise, a = 0.1. 
The existence of noise increases the correlation in neural firing at non-zero 
latencies, compared to the noise- free case (dashed line). To understand this 
seemingly counter-intuitive result, note that the spike correlation function 
in the noise-free case is identical to the spike auto-correlation function of a 
single neuron, and includes a delta-function contribution at At = 0. When 
the two neurons receive independent noise, they no longer fire precisely a 
the same time, leading to a broadening of the delta function into a peak of 
finite width. 

Panel b of the figure shows the spike correlation function of two neurons 
with different thresholds, 9\ = and 02 = a, for several values of a: (no 
noise), 0.4, and 1. Increasing a reduces the correlation between the two 
neurons. Note, however, that the preference of neuron 2 to fire after neuron 
1 is clearly evident at time scales comparable to r, even when the noise and 
signal have the same standard deviation. 

The inset shows c(At) fr\ri — 1 in a case where the noise is much larger 
than the signal, a = 10. Here c(At) is very well approximated using a 
linearization with respect to the cross- correlation function (dotted line), as 
described in Appendix B, Eq. (63). This approximation works quite well 
even if a is of order unity, as can be seen in the main plot (Fig. 6b, dotted 
line.) 
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6 Comparison with integrate-and-fire neurons 



Having characterized the spike correlation statistics of the threshold crossing 
model, we may ask whether the results carry over to other models of neural 
firing. In even the most simple leaky integrate-and-fire (LIF) model an 
analytical form is not known for the spike correlation functions. Hence the 
discussion in this section is based on numerical simulation of model neurons, 
using several variations of the LIF model. We focus on the relative timing of 
spikes elicited by neurons that receive the same stimulus but differ in their 
thresholds. 

In the following examples, where we consider pairs of neurons that differ 
in 6, we suppose that the neurons are identical, but differ in the mean value 
of their fluctuating input currents. Hence, after shifting the membrane po- 
tential to compensate for the different baseline currents, the reset potentials 
u r in the two neurons differ, but the difference 6 — u r is fixed. 

Figure 7 a shows the spike correlation function of two such neurons, 
evaluated from simulation of model neurons according to the LIF model 
of Eqs. (25)-(26), with 6>i = 0.8a, 6 2 = a, and 9 - u r = 0.5a. The only 
difference between the LIF model and the threshold-crossing model is the 
existence of a reset in the membrane potential following each spike. Nev- 
ertheless, comparing Fig. 7 a with Fig. 4 a reveals that there is very little 
resemblance between the spike correlation functions obtained from the two 
models. Most notably, in the LIF model there is almost no preference for 
one of the neurons to fire later than the other one. 

To better understand this discrepancy, Fig. 7 b shows an example of spike 
trains generated by the two neurons in the LIF model (top two traces), and 
in the threshold crossing model (bottom two traces). Compared to the LIF 
model, which tends to produce bursts of spikes, the threshold crossing model 
tends to generate much more isolated spikes. 

The inset in Fig. 7 b shows an example of bursts of spikes, generated 
by the LIF neurons after a relatively long silent period (top two traces). 
The first spike generated by the neuron with the smaller threshold (black) 
precedes the first spike generated by the neuron with the higher threshold 
(gray), preserving the tendency that is observed in the threshold crossing 
model. This tendency is typically observed also in other isolated bursts. 

Within the bursts, however, there is no clear relative timing of spikes 
of the two LIF neurons, because a spike in one neuron can be paired with 
a spike in the other neuron that either precedes it or comes after it. This 
suggests that the existence of bursts masks the tendency of the neurons 
to fire in a particular order. Motivated by these observations, we consider 
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several situations in which the LIF model can exhibit a characteristic order of 
spike timing, in similarity to the prediction of the threshold-crossing model. 

Sparse firing 

In situations where the LIF neurons fire isolated spikes, we may expect 
a clear order of spike timing to emerge, and to be reflected in the spike 
correlation function. 

Sparse firing due to a large potential reset, compared to the input vari- 
ance. With increase of 9—u r , compared to <r, spikes tend to become more iso- 
lated. With fixed biophysical parameters of the cell, an increase in (9—u r )/a 
corresponds to a decrease in the standard deviation of the fluctuating in- 
put current (see, also, Appendix D). Figure 8a shows that, as (0 — u r )/a 
increases in both neurons, a pronounced asymmetry develops in their spike 
cross-correlation function. 

Sparse firing due to refractoriness. Sparse firing may result, alterna- 
tively, from the existence of an additional refractory mechanism. There 
are many possible ways to model refractoriness in a LIF model (Gerstner 
&; Kistler, 2002), and here we consider one such possibility. After each 
spike, the membrane potential u is kept fixed at u r during a waiting period 
whose length varies randomly from spike to spike. After the waiting pe- 
riod the membrane potential continues to evolve according to Eq. (25), and 
the distribution function of waiting times decays exponentially with a time 
constant r r . Spike correlation functions, obtained from simulations of this 
model with r r = 2r, are shown in Fig. 8 b for several pairs of thresholds, 
exhibiting a clear asymmetry in firing order. Note that in this example, the 
membrane time constant is shorter than the typical refractory time r r . 

Bursting events 

Neurons often generate distinct bursts of spikes, sometimes referred to as 
events, when presented with stimuli possessing rapid temporal fluctuations 
(Berry et al., 1997; Berry & Meister, 1998b, 1998a). It has been argued that 
the timing of events, as marked by their first spike, can convey significant 
information about the stimulus (Hopfield, 1995; Gollisch & Meister, 2008). 
Hence it is of interest to understand the factors influencing the relative 
timing of these spikes. In the following, we consider a model where bursting 
events are more clearly separated from each other than in the simple LIF 
model. We then interpret the predictions of the threshold-crossing model in 
relation to the timing of burst onsets. 
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To reproduce well-separated spiking bursts in a neural model, we con- 
sider a generalization of the LIF model, introduced by Keat et al. (2001). 
This model was shown to reliably predict the structure of spike trains gen- 
erated by retinal ganglion cells from several different species, in response to 
rapidly flickering stimuli. In a version of this model without noise, a gener- 
ating potential g(t) is related to the stimulus by Eq. (5), as in the threshold 
crossing model. Spiking occurs when g{t) crosses a time-dependent thresh- 
old b{t) that increases following each spike, 

b{t) =0 + B f x(Oexp i^~f) ( 41 ) 

We refer to this model as the variable-threshold (VT) model. 

We may interpret g — b + as the neuron's membrane potential: with 
this interpretation the discrete increase in b following a spike corresponds 
to a reset of the membrane potential from to 6 — B. The VT model is 
precisely equivalent to the LIF model of equations (25)-(26) if / is a double 
exponential filter, as in Eq. (8), and if t v = t\. For simplicity we take / to 
be an 'alpha' filter, Eq. (11). However, we choose the time scale of threshold 
recovery t p to be larger than the neuron's membrane time constant, r p = 5r. 
This ratio roughly matches the relation between t p and the shape of filters 
that were found in Keat et al. (2001) to provide a good description of spike 
trains from retinal ganglion cells. 

A spike train generated by this model is shown in Fig. 9b, and consists 
of clearly separated events. The red lines in Fig. 9 b represent spikes gener- 
ated by a threshold-crossing model with the same generating potential and 
threshold as in the VT model. These spikes roughly match the onset of bust- 
ing events in the VT model. The threshold-crossing model thus provides a 
coarse-grained description of events generated by the VT model. 

We next isolate the first spikes of events in the VT model by discard- 
ing any spike that occurred within a time delay D = 2r from a previous 
spike 7 . Fig 9 c shows the cross-correlation function of these specifically 
selected spikes in several neurons with thresholds 0/a = 0.5, 0.2, —0.2, 
and —0.5, that were all presented with the same stimulus. These cross- 
correlation functions display a strong asymmetry between positive and neg- 
ative At, and are qualitatively very similar to the spike cross correlations in 
the threshold crossing model, shown in Fig. 5. In contrast the spike cross 
correlation functions calculated directly from all the spikes, Fig. 9 a, are very 
different from the prediction of the threshold-crossing model. 

7 Because the events are clearly separated clearly from each other, the results are in- 
sensitive to the precise choice of D. 
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These results suggest that the onsets of bursting events in neurons that 
receive the same stimulus, but differ in their threshold, can exhibit strong 
asymmetry in their relative timing, and that the correlation functions of 
these event onsets can be qualitatively described by the simplified threshold- 
crossing model. 

7 Summary 

We considered in this work a relatively simple deterministic process that pro- 
duces a discrete spike train from an analog, continuous signal. The timing of 
spikes in this model is more precisely controlled by the stimulus than typi- 
cally predicted by rate-based models. This aspect of the model is motivated 
by observations of precisely timed spiking in neural assays, particularly in 
response to stimuli that possess strong temporal modulations (Mainen & 
Sejnowski, 1995; Berry et al., 1997; R. R. de Ruyter van Steveninck et al., 
1997; Meister & Berry, 1999; Uzzell & Chichilnisky, 2004). Some of the 
salient properties of spike correlation functions in this model are summa- 
rized below. 

Most notably, two neurons receiving identical inputs can show a pref- 
erence for one neuron to fire later than the other, although there is no 
monosynaptic connectivity between them (Figs. 4a and 5.) This preference 
is prominent even if the thresholds of the two neurons are similar (but not 
equal), and if there is only partial overlap in their input (Fig. 6.) 

Comparison with the leaky integrate-and-fire model shows that bursting 
in LIF neurons often masks the preference of neurons to fire in a particu- 
lar order (Fig. 7). However, there are several situations in which such a 
preference may be observed in LIF neurons. First, LIF neurons can dis- 
play a preference to fire in a particular order if they produce sparse spiking, 
e.g., due to a large hyperpolarizing step in the membrane potential after 
each spike, compared to the standard deviation of the stimulus, or due to 
other sources of refractoriness (Fig. 8). Second, when the neurons generate 
clearly separated bursting events, the first spikes of these events have cross- 
correlation functions that are similar to those predicted by the threshold 
crossing model (Fig. 9). 

The spike cross-correlation function of neurons that differ only in their 
thresholds may be probed experimentally by repeatedly presenting the same 
stimulus to a single neuron, while injecting varying amounts of current into 
the neuron from trial to trial. Correlations between different trials can then 
effectively measure the spike correlation function of two identical neurons 
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with different thresholds. This approach was recently taken in (Markowitz, 
Collman, Brody, & Hopfield, 2008), where a Gaussian stimulus, mimicking 
the spectral properties of Gamma oscillations, was injected into rat pyrami- 
dal neurons from the somatosensory cortex. While spike cross-correlation 
functions are not shown in (Markowitz et al., 2008), spike trains from differ- 
ent trials exhibit a clear modulation of spike timing by the varying injected 
current, suggesting that a strong asymmetry exists in the spike-correlation 
function of neurons with different thresholds. Furthermore, the characteris- 
tic latency between the firing of two neurons appears to increase linearly as 
a function of the difference in injected currents. This roughly linear depen- 
dence persists over a wide range of current differences, in similarity to our 
results from the threshold-crossing model. 

We also considered in this work the spike-triggered average stimulus 
(STA). The STA in our model is a linear combination of the filter and its 
derivative. A simple geometrical argument shows that this result extends 
to a larger class of models: it should hold whenever the stimulus is passed 
through a linear kernel, if spike decisions are then based strictly on the 
output of the kernel and its derivative. Varying the threshold modifies the 
relative weight of the filter and its derivative (Fig 1). It will be interesting 
to probe for such a dependence of the STA on threshold in real neurons. 
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Note added in proof 

Tchumatchenko et al. have recently considered a threshold-crossing model, 
similar to the one presented in this work. A preprint of their work has been 
made available on the arXiv.org e-Print archive while our manuscript was 
in review (Tchumatchenko, Malyshev, Geisel, Volgushev, & Wolf, 2008). 
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A Spike triggered average stimulus and covariance 



To avoid infinities in the calculation we assume first that s is a correlated 
Gaussian process with a correlation function w s (At). At the end of the 
calculation we take the limit 



w s (At) -> ald{t - t') 
The correlation function of g is 

w(At) = J dt J dt' f(t)w s (t -t' + At)f(t') 



(42) 



(43) 



To evaluate /sta we need the joint probability distribution function for 
g(0) = 9, g(0) = q±, and s(At) = s 2 , in terms of which /sta is given by 



roc poo 

/sTA(At) = - / ds 2 / dqis 2 -qi-p(0,q 1 ;s 2 ) 

r J-oo JO 

The joint probability distribution of p(6, q±; s 2 ) is given by 



P(0,qi;s 2 ) = -^exp 



2 



(44) 



(45) 



where Q T = (9, qi ,s 2 ), Z = [(2vr) n deU] 1/2 , and 



w(0) (fow s )(-At) 

-w"(0) (f'ow s )(-At) 

(fow s )(-At) (fow s )(-At) w s (0) 



(46) 



Evaluating the integral (44) yields 

/ STA (At) = ^(/o«, s )(-At) + 



7T 



-2w"(0) 



(f'o Ws )(-At) (47) 



In the limit where w s is uncorrelated, Eq. (42), we get Eq. (24). 
The spike-triggered covariance C(t\,t 2 ) is defined as 

C(h,t 2 ) = l - {[ 8 (h) - /sta(*i)] [s(h) - /sta(* 2 )] • X(0)> 
A calculation similar to that outlined for the STA yields 



(48) 
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The eigenfunctions of the covariance operator and their corresponding 
eigenvalues Aj are determined by the equation 



/ 



C(t 1 ,t 2 )Mt2)dt 2 = XiMh) (50) 



From Eq. (49) we see that there are only two eigenfunctions ^2 (At) with 
eigenvalues that differ from Oq. (The significance of these eigenfunctions 
is that the variance of the stimulus's projection on ipip, when conditioned 
on the occurrence of a spike at t = 0, is different from its nominal value 
of 00.) The first eigenfunction tpi(At) is proportional to /(—At). Because 
the generating potential must be equal to the threshold at t = 0, when a 
spike is produced, the variance of the stimulus's projection on / must be 
zero: Indeed, Ai = 0. The second eigenfunction i[) 2 (At) oc /'(—At), and 
A 2 = (2 - vr/2)a2. 



B Spike correlations 



We consider two model neurons that may differ in their generating potentials 
<7i )2 and in their thresholds 0\ t2 . The generating potentials g± t 2 are assumed 
to be jointly Gaussian and stationary, characterized by correlation functions: 

Wij (At) = { 9i (0) gj (At)) (51) 

The stationary nature of g± j2 implies that wu(— At) = wu(At) and that, 
similarly, W22(—At) = W22{At). For the cross-correlation functions it only 
implies that w\ 2 (—At) = w 2 \(At). In the rest of this appendix we assume 
that, in addition, 

^12 (At) = ^21 (At) (52) 

which is correct throughout Sec. V. In this case the expressions for the 
joint probability distribution and for the spike correlation function simplify 
considerably. 

We first need to evaluate the joint probability distribution for g\ (0) = 0±, 
9i (0) = 9i, 92 (At) = 2 , and 52 (At) = q 2 . This is given by 



p{Oi,qv,02;q2) = ^exp 



2 



(53) 



where £ 5 



A 



(#1,91,02,92), Z 

w n (0) 


W12 (At) -w' 12 (At) 

w' 12 (At) -w'{ 2 {At) 



= [(2vr) n det^] 1/2 , and 

w 12 (At) 
-<i(0) -«/i 2 (At) 



1^22(0) 




w' 12 (At) 
-w'{ 2 (At) 


-^22(0) 



(54) 
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Since, in calculating the spike correlation function, 9\ and 9 2 are kept 
fixed, it is useful to re-express this quantity as a quadratic function of q\ 
and q 2 alone. We know that 



dqi j dq 2 p(9 1 ,q 1 ;6 2 ,q 2 ) = p(9i;6 2 ) 



(55) 



where p{9\]9 2 ) is the joint probability distribution of g\ (0) = 9±, g 2 (At) 

02, 



p(0i;0 2 ) 



i / ^ii(O)0f+U»2 2 (O)^-2^i2(At)0 1 02 

exp ' 



2t:D 



2D 2 



and 



£> = [w u (0)w 22 (0) -wj 2 (At)] 1/2 . 
Hence we can rewrite 51; 2 ,q 2 ) as 



(56) 
(57) 



(detm 1 ) 1 ^ 2 
P = p(0i;0 2 ) x ^ — -!■ exp 



2tt 



(q-qo) T "i ^q-qo) 



The matrix m and qo are found by collecting the quadratic and linear 
terms in qi, q 2 in Eq. (53), 



m 



and 



-<(0) -< 2 (At) 
-< 2 (At) -u&(0) 



[«/ 12 (At)]' 



w n {0) w 12 (At) 
wu(0)w 22 (0) - ^ 2 2 (At) L w 12 (At) w 22 (0) 

(58) 



qo 



w[ 2 (At) 



w 12 (At)d! - w u (0)6 2 
10^(0)0! - w 12 (At)e 2 



w u (0)w 2 2(0) -w\ 2 {At) 
Finally, we can write the spike correlation function as 



c(At)=p(9 1 ;9 2 )-I 



(59) 



(60) 



where 



(detm" 1 ) 172 [°° [°° f 1 xT -1 / 

J = — / dgi / dga gi ■ ga ■ exp --(q-qo) m (q - qo) 

(61) 



22 



Nearly independent inputs 



When the cross-correlation w\2 (At) is small compared to the auto-correlations 
u>n(0), 1^22(0), I and p can be expanded in powers of wuiA-t) and its 
derivatives. If wu^At) = the matrix m is diagonal, p(9i, q±; 62, 92) = 
p(@iiQi) ' p(@2,Q2), and c(At) = r\T2- To first order in w\2^At\ 



4 m 7 2v^| ^22(0) 



1/2 



[-^22(0)] 1/2 

wn(0) 



2 ^ w' 12 (At)+- 



(62) 

' ; ! / 1/2 

where Jo = [^11 (Oju/^O)] /(2tt)- By similarly expanding p(#i; #2) we find 
that 



c(At) — r±r 2 



6162 



nr 2 



+ 



^n(0)u;22(0) 



^12 (At) -X' 2 (At) 

2K' 1 (0K(0)] 1/2 



2 U-^' 2 (o)] i/2 u;22(o) [-w'^m^mm 



Identical inputs 



Equation (61) can be further simplified in the case where g\ and (72 are 
identical. In this case w\\ = w\2 = W21 = ^22 = w. The spike correlation 
function involves five parameters: w and its second derivative at zero, w, 
its derivative, and its second derivative at At. In Eq. (61), m and qo, 
Eqs. (58)-(59), are then 



m = 



-w"(0) -w"{At) 
-w"(At) -w"(0) 



[w'(At)Y 



w 2 (0) - w 2 (At) 



w(0) w(At) 
w(At) w(0) 



qo 



w'(At) 



w(At)0i - w(O)0 2 



w 2 {0) - w 2 {At) [ w{O)0i - w{At)0 2 J ' 



, (64) 
(65) 



and 

p(0i;0 2 ) 



2vr [w(0) 2 - w(At) 



. w(<d)(0l + 2 2 ) - 2w(At)0i0 2 
^ 2 ' ' XP 1 ~ 2 [^(0) 2 - u>(At) 2 ] 

(66) 

Note that from symmetry under time reversal we must have 



p(0i,qi;02,q2) = p(02, ~qi\ Oi, -qi) (67) 
This symmetry is reflected in the fact that (qo)i,2(#i, #2) = — (qo) 2 ,i(#2, 



23 



Because mn = rri22, the eigenvectors of m are 



1 



V2 V ±1 

and the corresponding eigenvalues are 



m± = — u/'(0) =F w"(At) 



[w'(At)Y 



(68) 



[-w(O) T w(At)] (69) 



w 2 (0) -w 2 (At) 

It is useful to rewrite / in the coordinates in which m is diagonal 

1 



u± = -j=(qi±q 2 ) 



in terms of which 



I = - I du+ / du_ (u 2 ^ — u 2 __) 
2 Jo 



(70) 



(71) 



27rmV 2 m 1//2 



exp 



(u + - u ,+) 2 (u_-u _) 5 



where 



1 



uo,± 



2m i 



w'(At) 



2m_ 



'(02TO!) 



^2w(0)Tvj(At) y " A ' " L/ (72) 

The inner integral can be expressed using the Gauss error function, yielding 
an expression for / that involves a single integral, 

(u+ - u ,+) 21 



r 

Jo 



du-j 



x < m_ (u + — uo,-)exp 
+m^ 2 (u + + iio,-)exp 



1 

47rm. 



1/2 

+ 



exp 



2m + 
(u+ +n ,-) 2 

(u + - u _) 2 



(73) 



+\/o ( m - + u o- -«+) 



erf 



2m_ 

V / 277jT 



erf 



u _ + u + 



The limit At -> oo 

Assuming that u> and all its derivatives at t = At tend to zero when At — > oo, 
the probability distributions at t = and at t = At decouple, and to leading 
order c(At) — > r\V2- 
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To evaluate the deviation from independent spiking, we can use equation 
(63) in the case where w±i = W22 = ^12 = ^21 = w: 



c(At) — r\T2 
r\T2 



+ 



6162 

Mo) 
1 

w{0) [-2w"(0)_ 



IT 

2w"(0) ' 

1/2 

(0 1 - 2 )w'{At) 



7T 



(74) 



The third term in this equation is antisymmetric in At and in 0\ — 62. 
Because w'(At) is typically negative for positive At, this term represents a 
small preference for the neuron with the higher threshold to spike after the 
neuron with the lower threshold. 



C The limit At 

We consider here the limit At — ► in the case of identical inputs. In this 
limit the matrix m becomes singular, requiring particular analysis in order 
to evaluate the correlation matrix. 

The argument of the exponential in Eq. (28) is maximal when q = qo 
[Eq. (65)] which approaches, when At — > 0, 

(S0)l, 2 - ^ (75) 

This is the derivative of g(t) if it follows a linear trajectory from g{0) = 0\ 
to g{At) = 02- However, this maximum is within the integration range in 
Eq. (28) only if {02—9i)/At is positive, i.e., only if the neuron with the higher 
threshold spikes after the neuron with the lower threshold. Accordingly, the 
behavior of c(At) is substantially different for positive and negative At. 

For simplicity of the notation in this appendix, we assume that At > 
and that 02 — Q\ may be either positive, negative, or zero. To relate this to 
the presentation in Sec. V, where we assumed that 02 — 0\ is positive, recall 
that c(— At) is the same as c(At) with the two thresholds 0\,02 exchanged. 

To evaluate the behavior of the integral / at small At we need to expand 
the expressions for m+, m_, u+, and u_ in this limit, where we use the 
representation of / in Eqs. (71)-(72). The leading terms in these expansions 
are shown in Table 1 in two cases: (i) W3 / 0. (ii) All W{ with odd index 
vanish. 
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VF 3 = W 5 = 










W 3 At 








V2u_ 







Table 1: Leading order terms in the expansion of m± and u± for small At. 
If 61 = 62, u + = to all orders. 



The case 2 - Q\ > 



This case is most simply treated in the original (q) coordinates. As At — > 
(90)1,2 ~ (At)" 1 whereas the standard deviation of the Gaussian in Eq. (28) 
tends to zero as y/At in both of the principal directions (or faster if W3 = 
0). We can therefore replace the lower integration limits by — 00 with an 
exponentially small error. Furthermore, the integral in Eq. (28) becomes 
dominated by the peak of the Gaussian, as it becomes sharper with the 
decrease of At. 

To leading order in At the integral is thus simply 



(<7o)i • (90)2 



At 2 



To evaluate p(9±; 62) we note that, in Eq. (66), 

w 2 {0) - w 2 (At) ~ -W W 2 (At) 2 
and by expanding the argument of the exponential we obtain 



p(0i;O 2 ) 

where 



27r(-WW 2 ) 1 / 2 |At| 



x exp 



1 f9o-6i 



2W 2 



At 



+ 



7i_ 
At 



To 



7o 



72 Wf 



(AW 2 



3W 2 W 4 )(9 2 - 9i)' 



8W 



(9 1 + 9 2 ) 2 . 



(76) 



(77) 



(78) 



(79) 
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and 



7i 



W 3 
'6W 2 2 



(02 - Oxf 



(80) 



The leading contribution in Eq. (78) is proportional to the probability den- 
sity to have a derivative equal to (62 — 9\)/At. Combining Eqs. (76) and 
(78) we obtain the small At behavior of c(At), Eq. (32). 

When 62 — 6\ — > the peak in the spike correlation function gradually 
becomes the delta- function contribution to the spike auto-correlation func- 
tion, Eq. (29): For sufficiently small 6 2 - 9i, and for At < -W 2 /W 3 , the 
spike correlation function is approximately 



c(At) 



(0 2 -9 x f 



2Tr(-W W 2 ) 1 / 2 At :i 



x exp 



1 



2W 2W 2 V At 



6o-9 



(81) 



The integral of this function from to infinity is equal to r, whereas the 
width of the peak and its position both tend to zero when 92 — 9\ — > 0. 
Hence the small At behavior of the spike correlation function approaches, 
when 02 — 9i — > 0, 

rS(At). (82) 



The case 9 2 = 61 

The maximum of the Gaussian in Eq. (71) is at u + = and U- ~ 0At 
(where 9 = 9\ = 9 2l ) which lies outside the range of integration [The (+,+) 
quadrant in the q coordinates.] However, if W3 / the distance from the 
(+,+) quadrant is small compared to the standard deviation of the Gaussian 
(~ \/At). As a result, the integral in (71) can be treated as if uq : ± = 0. 
The situation is different if W3 = 0, in which case the distance of q from the 
origin and the standard deviation in the «_ direction are of the same order 
of magnitude. 

(i) W3 / 0. In this case we have, to leading order in At, 

0,2 r.,2 



I ~ 
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(83) 



(84) 
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Hence the spike correlation function tends to a finite value when At — > 0, 
Eq. (30). 

(«) W3 = W5 = 0. To treat this case we rescale u± by the leading depen- 
dence on At in the expansion of (m±) 1 l 2 (this procedure can also be applied 
in the case W3 / 0, in order to derive Eq. (83) in a more formal manner). 
We use the following notation, 



m + ~ a + (At) 4 
u + = u+ ■ (At) 2 



m,- ~ a_(At) 2 
u_ = u_ • At 



u_ ^ «At (85) 



where 



a = — 



72 \ W 2 
W 2 



W 6 
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+ W 4 
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(86) 




du_ [(At) 4 ^ - (At) V 



(87) 



Using Eq. (84) we see that, in contrast to the case W3 7^ 0, the spike 
correlation function tends to zero when At — > 0, 



c(At) ~ (At)' 



(88) 



The case 9 2 - 9 1 < 

In this case, as in the case 62 > #1, (<?o)i,2 ~ At -1 , but now qo is outside 
the integration region, far from it on the scale of the Gaussian's standard 
deviation. We denote 



u+,o ^ - 



a 
At 



(89) 



where a = \f2(0\ — 62) > 0. It is then seen from Eq. (71) that / decays to 
zero exponentially, with a leading contribution of the form 



log(J) ~ - 



a 



2m+{At) 2 
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m + 



0\ — 62 
At 



(90) 
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which corresponds, if W3 ^ to 



and, if W 3 = W 5 = 0, to 



(e 1 -e 2 ) 2 



lo g (/)^-^4f + - (92) 

In contrast to / in the case 62 — 6\ > 0, which diverges when At — ► 0, in 
the case 62 — 6\ < / strongly decays to zero. A more precise treatment of 
this integral requires evaluation of additional terms in the expansion shown 
in Table 1. 

The probability p(0±; 62) is the same as in the case 62 — 61 > 0. Like /, p 
decays exponentially when At — > 0, but the leading power of (At) -1 inside 
the exponential is smaller, so it becomes relevant only if log/ is expanded 
beyond the leading order. 



D Parameterization of the LIF model 

The notation in equations (25)-(26) was chosen for mathematical simplicity. 
To see how the parameters in these equations are related to biophysical prop- 
erties of the neuron, we start with a standard equation for the membrane 
potential dynamics of a LIF neuron (Dayan & Abbott, 2001), 

dv 

n— = E - v + R m I c (93) 

where Eq is the resting potential of the cell, R m is the membrane resistance, 
and I c is the input current. Further, the potential v resets to v r whenever 
it reaches a threshold vq. 

To reparametrize this equation as in Eq. (25), we first define 

/ = R m (I c - </ c » (94) 

where {I c ) is the mean (baseline) value of the fluctuating current I c . Due to 
the multiplication by R m , I has the same units as the potential u. Second, 
we shift the potential v by defining: 

u = v-E -R m (I c ) (95) 
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With these definitions we obtain Eq. (25), repeated here, 

n^ = -u + l (96) 

where / has zero mean, and 9 — u r = vg — v r . 

The properties of the spike train depend on the ratio between 9 — u r and 
the standard deviation a. The following example demonstrates what is a 
reasonable range of values for this ratio. The standard deviation of / is equal 
to \[2o (assuming that T\ = r-i). Suppose that the standard deviation of the 
current I c is 0.1 nA. Assuming that R m = 100 Mil (which corresponds to a 
neuron with surface area 10 -4 cm 2 and membrane time constant t\ = 10 ms), 
and that the difference between the threshold membrane potential and the 
reset potential is 10 mV we get, using Eq. (94), 9 — u r ~ 0.7a. 

The ratio (9 — u r ) jo can vary considerably depending on the biophysical 
parameters of the cell and, more importantly, depending on the standard 
deviation of the fluctuating current. A smaller standard deviation of the 
fluctuating current corresponds to a larger value of (9 — u r )/a. 



E Correlations in the LN model 

In the LN model the spike times are an inhomogeneous Poisson process, 
with a rate given by 

r(t) = <*[(/ <>*)(*)] (97) 

where eft may be regarded as the transfer function of the neuron, if we think 
of s as a direct input to the neuron (more generally, we may think of Eq. (97) 
as a phenomenological relation between stimulus and firing rate.) 

If there was no nonlinearity [<p(x) = x] we would have (r(t)r(t')) ~ 
w(t — t'), where w is the auto-correlation function of / o s. More generally, 
we need to take into account the particular form of <fi. We focus here on the 
case where two neurons, labeled 1 and 2, receive the same input and have the 
same filter, but possibly differ in the nonlinear function <f>. As in previous 
sections we assume that the stimulus s is Gaussian (and so is g = fos). We 
then have, 

/OO /'OO 
dgi / dg 2 M9i)M92)p(9i;92) (98) 
-oo J —oo 

where p(gi;g2) is given by Eq. (66). The spike correlation function thus 
depends on w(0) and on w(At): in comparison, in the threshold crossing 
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model it depends also on the first and second derivatives of w at these two 
points. 

For example, if 

<hM = A> { ° * ' a > I 1 " <") 

[ fl 1 - #1,2 , 9 > Vl,2 

the spike correlation function is 

roc roc 

c(At) = $ / dsi / d 52 ( ffl - ei)(g2 - 2 )p{gi,92) (100) 

Note that the spike correlation function in the LN model is symmetric in 
At — > —At despite the different offsets for the two neurons. 

In the special case where 9\ = 62 = 0, the spike (auto-)correlation func- 
tion can be written in a relatively simple form, 

. . . 1 r^r, 7TT-r-r w(At) \w(At)\ \w(At)\ 

c(At) = y / w ^0)- w ^At) + ^-^ + LA^ arctan I ^ 

27T 4 2lT y/w 2 (0) ~ ^ (At) 

(101) 

For any two thresholds, the behavior for large |At| can be found from 
the expansion of p(<7i; 32)) from which we obtain 



c(At) ~ r iro 2 + ^^(At)erfc f / ° 1 ) erfc ( / 2 ) 



(102) 



where 

r , = (- 4-) - ^erfc f * ) (103) 

are the firing rates of the two neurons and erfc is the complementary error 
function. 

For a more direct comparison with the threshold crossing model it is 
instructive to consider a transfer function that is concentrated at a particular 
value 9: 4>(x) = 4>q5{x — 9). We then find that for large |At|, 

~ W ° - (104) 

which is the same as the first term in Eq. (74). The firing rate in the LN 
model, appearing in the left hand side of Eq. (104), is given by 



r0,i y/^^O) 6 ^ 



2iu(0) 



(105) 
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Figure 1: (a) STA of a model neuron spiking in response to upward threshold 
crossing of the generating potential described by Eqs. (5), (11), and (12). 
The STA, Eq. (24), is plotted for three values of the threshold: 9 = 2.5a, 
a, and 0, where a is the standard deviation of g(t) (solid line: analytic 
expression, dotted lines: simulation), (b) A comparison of the STA for 
6 = 2.5a (solid line) with simulation results from a leaky integrate-and-fire 
neuron with time constants t\ = T2 = r and the same threshold, receiving 
s as its input, and resetting its membrane potential to after each spike 
(dotted line). The dashed line shows the first term only of Eq. (24). (c) A 
similar comparison as in (b), with 6 = 0.25(7. (solid line, threshold crossing 
model, dashed line, leaky integrate-and-fire model.) 
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Figure 2: Spike autocorrelation function, c(At), divided by the average 
firing rate r, for three values of the threshold: 9 = 0, 2a, and 3o\ The 
generating potential is described by Eqs. (5), (11), and (12). (b) A simi- 
lar plot for a generating potential having a Gaussian correlation function, 
w(t) = a 2 exp(-t 2 /2r 2 ). 
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a Af > b At < 




Figure 3: Spike generation events in the threshold crossing model, for two 
neurons with thresholds 0\ (thick circle) and 02 > 9± (thin circle). In the 
left plot neuron 2 fires after neuron 1, whereas in the right plot neuronl 
fires after neuron 2. In the latter case the generating potential must reverse 
the sign of its derivative twice within the time interval separating the two 
spikes. 
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Figure 4: (a) The spike correlation function, c(At)/r±r2 (solid line,) of two 
model neurons with thresholds 9\ = 0.8a and 62 = <J, firing in response 
to the same stimulus (same as in Fig. 2a). The dashed line shows the 
approximation for small At, Eqs. (76) and (78)-(80). At large |At| the two 
neurons become decorrelated, and c(Ai)/rir2 approaches unity. The inset 
shows the approximation for weakly correlated neurons, Eq. (74) (dashed 
line) compared with the actual correlation function (solid line), (b) The 
prediction of a LN model for the spike correlation function, c(At) /r±r2, 
shown for two different forms of the non-linearity, as described in the text 
(linear rectification, solid line; delta function, dashed line). 
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Figure 5: Spike correlation function, c(Ai)/nr2, for two model neurons with 
02 = 0.5o" and 0\ = 0.5(7 (dashed line), 0.2(7, — 0.2a, and — 0.5(7 (solid lines). 
The generating potential is described by Eqs. (5), (8), and (8). Arrows point 
to the position of the peak according to the prediction of Eq. (34). 
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Figure 6: (a) Spike correlation function, c(At) /r±r2, of two model neurons 
with identical thresholds 6 = a, receiving an identical stimulus and uncor- 
rected noise, Eq. (38) (solid line.) The standard deviation of the noise, 
divided by that of the common stimulus, is a = 0.1 (other parameters are 
described in the text.) This is compared with the case where there is no 
noise, a = (dashed line.) (b) Spike correlation function of two neurons 
receiving similar input as in (a), but differing in threshold: 6\ = and 
62 = <J, where a is the standard deviation of the common stimulus s. The 
correlation function is plotted for three different ratios between the noise 
standard deviation and the stimulus's standard deviation: a = 0, 0.4, and 
1. The inset shows c(At)/r\T2 — 1 in a case where the common stimulus is 
weak compared to the noise, a = 10 (solid line). The dotted line in the inset 
and in the main plot (for a = 1) is the approximation for the case of weak 
correlation, Eq. (63). 



40 



At/x 



■mi i in 



I I 



i i i 



1 







— n 


lillllllll 1 








46 


48 

t/x 



20 40 60 

t/x 



Figure 7: (a) Spike cross-correlation function of two simulated LIF neurons, 
Eqs. (25)-(26) with n = T2 = r, and with thresholds 0.8 a and 1.0 a, receiv- 
ing an identical stimulus (as described in the text, up to a baseline shift). 
In both neurons 6 — u r = 0.5a. (b) Spike trains generated by the two LIF 
neurons (top two traces: black, 9\ = 0.8 and gray, 62 = 1) and by thresh- 
old crossing model neurons responding to the same stimulus (bottom two 
traces). Inset: blow-up of an interval showing a burst generated by the two 
LIF neurons. 
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Figure 8: (a) Spike cross correlation functions of LIF neurons with increasing 
values 6 — u r , the hyperpolarizing jump in potential following each spike 
(alternating solid and dashed lines: 0.5a, 2a, 4a, and 6a. In all traces the 
thresholds are 9\ = 0.8 and gray, 62 = 1. The trace for 6 — u r = 0.5<r is 
identical to the one in Fig. 7 (a), (b) Spike cross-correlation functions of 
pairs of LIF neurons with refractoriness (as described in the main text.) 
In all pairs the second neuron has a threshold 62 = 1-0, whereas the first 
neuron has a threshold 9\ = 0.2 (solid line), 0.5 (dashed line), and 0.8 (solid 
line), and 6 — u r = 0.5a. 
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Figure 9: (a) Spike cross-correlation in the model described in Sec. VI. The 
generating potential is described by Eqs. (5), (8), and (9). The thresholds 
are 9\ = 0.8cr and 62 = a, t p = 5r and B = 0.1. (b) An example of a 
spike train generated by this model (short vertical lines). The threshold 
6 = a. For comparison, the tall vertical lines represent spike times in the 
threshold- crossing model. The generating potential, common to both mod- 
els, is plotted in gray. Top right: blow-up of a shorter interval, showing 
three spiking events, (c) Cross correlation of event times in the same model, 
between a model neuron with threshold 9\ = 0.5 and three model neurons 
with thresholds 62 = 0.2, —0.2, and —0.5. Other parameters are as in pan- 
els a and b. Event onsets are isolated by discarding any spike that occurs 
within a time frame of 2r from a previous spike. 
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